/*
 * @Description: 优化法估计矩阵1范数
 * @Autor: kaikaima
 * @Date: 2021-04-06 20:27:16
 * @LastEditors: kaikaima
 * @LastEditTime: 2021-04-17 12:26:40
 */

#include <iostream>
#include <cmath>
#include "Linear_equations.h"

int main(int argc, char* argv[]){
    for(int N=5;N<=30;N++){
        double key1=0,key2=0,rnorm=0,bnorm=0,_xnorm=0,xnorm=0;
        double z[N],w[N],v[N];
        L_Equations P75_2(N);

        double randomx[N];
        for(int i=0;i<N;i++){
            std::cin>>randomx[i];
            if(xnorm<fabs(randomx[i])) xnorm=fabs(randomx[i]);
        }

        double* A[N+1];
        for(int i=0;i<=N;i++) A[i]=new double[N];

        for(int i=0;i<N;i++){
            A[N][i]=1.0/N;
            for(int j=0;j<N;j++){
                if(i==j) A[i][j]=1.0;
                else if(i>j) A[i][j]=-1.0;
                else if(j==N-1) A[i][j]=1.0;
                else A[i][j]=0;
            }
        }

        //求A的无穷范数
        for(int i=0;i<N;i++){
            v[i]=0;
            for(int j=0;j<N;j++)
                v[i]+=fabs(A[i][j]);
            if(key1<v[i]) key1=v[i];
        }

        //解方程得r
        {
            P75_2.set(A);

            double* _b=P75_2*randomx;
            for(int i=0;i<N;i++){
                if(bnorm<fabs(_b[i])) bnorm=fabs(_b[i]);
            }

            P75_2.set(_b);
            double* _x=P75_2.Gauss_solve();
            for(int i=0;i<N;i++)
                if(_xnorm<fabs(_x[i]-randomx[i])) _xnorm=fabs(_x[i]-randomx[i]);

            double _r[N];
            for(int i=0;i<N;i++){
                _r[i]=0;
                for(int j=0;j<N;j++)
                    _r[i]+=(A[i][j]*_x[j]);
                _r[i]=_b[i]-_r[i];
                if(rnorm<fabs(_r[i])) rnorm=fabs(_r[i]);
            }

            delete _b;
        }

        while(true){
            {
                P75_2.set(A);
                P75_2.transposition();
                double* _w=P75_2.Gauss_solve();
                for(int i=0;i<N;i++) w[i]=_w[i];
            }

            for(int i=0;i<N;i++)
                if(w[i]>0) v[i]=1;
                else if(w[i]==0) v[i]=0;
                else v[i]=-1;
            {
                P75_2.set(A);
                P75_2.set(v);
                double* _z=P75_2.Gauss_solve();
                for(int i=0;i<N;i++) z[i]=_z[i];
            }
            
            //get ||z||∞ zTx
            int zn=0;
            double zx=0;
            for(int i=1;i<N;i++) if(fabs(z[i])>fabs(z[zn])) zn=i;
            for(int i=0;i<N;i++) zx+=(z[i]*A[N][i]);

            if(fabs(z[zn])<=zx){
                for(int i=0;i<N;i++)
                key2+=fabs(w[i]);
                P75_2.set(A);
                P75_2.transposition();
                break;
            }
            else{
                for(int i=0;i<N;i++)
                    if(i==zn) A[N][i]=1;
                    else A[N][i]=0;
            }
        }

        std::cout<<N<<"阶矩阵的估计精度: "<<key1*key2*rnorm/bnorm<<std::endl;
        std::cout<<N<<"阶列主元运算实际精度: "<<_xnorm/xnorm<<std::endl<<std::endl;

        for(int i=0;i<=N;i++)
        delete A[i];
    }
    return 0;
}